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ABSTRACT 

We investigate asymptotic convergence in the Ax — > limit as a tool for determining 
whether numerical computations involving shocks are accurate. We use one- 
dimensional operator-split finite-difference schemes for hydrodynamics with a von 
Neumann artificial viscosity. An internal-energy scheme converges to demonstrably 
wrong solutions. We associate this failure with the presence of discontinuities in the 
limiting solution. Our extension of the Lax-Wendroff theorem guarantees that certain 
conservative, operator-split schemes converge to the correct continuum solution. For 
such a total-energy scheme applied to the formation of a single shock, convergence of a 
Cauchy error approaches the expected rate slowly. We relate this slowness to the effect 
of varying diffusion, due to varying linear artificial-viscous length, on small-amplitude 
waves. In an appendix we discuss the scaling of shock-transition regions with viscous 
lengths, and exhibit several difficulties for attempts to make extrapolations. 



1. Introduction 



Finite-difference hydrodynamics schemes find widespread use in astrophysical applications, 
where they are often parts of computations that involve other significant physical processes, 
such as external heating and radiative cooling. In this paper we consider numerical methods for 
establishing the accuracy of a calculation in the presence of shocks. Our primary tool is testing 
for convergence. We require that calculations not merely converge, but converge at or near 
the expected asymptotic rate. This accuracy criterion should readily lend itself to automated 
refinement methods in which grids are made denser only when and where required for accuracy. 

For our investigations we use the one-dimensional Euler equations modified by the inclusion 
of an explicit artificial viscosity. The numerical schemes are operator-split and explicit. We 
describe two schemes: a non- conservative form that uses the internal-energy density as a primary 
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variable, as is common in astrophysical calculations (e.g., Norman & Winkler 198f:, [Stone & 



Norman 1992, and [Hawley, Smarr &: Wilson 198'j| ), and a similar conservative form that instead 



uses the total-energy density as a primary variable. 

In this paper we investigate the performance of these schemes in simple test problems 
with shocks. After outlining the numerical methods in section 2, we show in section 3 that the 
non-conservative internal-energy scheme converges to incorrect solutions when shocks develop. 
This failure occurs when artificial-viscous lengths and gridsize are decreased proportionally, as is 
the usual practice. When viscous lengths are held fixed as the gridsize tends to zero, however, 
the scheme converges to a correct solution (albeit to a different physical problem). Section 4 
discusses the performance of both schemes in a shock-tube problem. We show explicitly that the 
internal-energy scheme (for a particular choice of viscous lengths) converges toward a solution with 
shock speed too low by approximately 0.2%. The other, conservative difference scheme converges 
toward the correct solution even in the presence of the shock. (We show in Appendix A that its 
correct behavior is a consequence of the Lax-Wendroff theorem [[Lax &: Wendroff 1960 1, which we 
extend to operator-split schemes.) 

Then in section 5 we examine rates of convergence for the conservative scheme. We show that 
the presence of linear artificial viscosity (proportional to dv/dx) introduces a substantial amount 
of diffusion. In sequences in which the artificial viscosity vanishes with the gridsize, the variation 
in diffusion has an observable effect on convergence rates for small-amplitude waves. Finally, we 
consider a wall-shock problem, which also contains such small-amplitude waves. We consider two 
kinds of sequences: (1) sequences in which the artificial viscosity enters as a constant term, and 
(2) sequences in which it varies with the gridsize. Holding the artificial-viscous lengths fixed leads 
quite readily to small changes and the expected convergence rate. In contrast, in the sequences 
where the viscosity vanishes with the gridsize, we observe that the convergence rate tends only 
slowly toward the rate expected. The effect of linear artificial viscosity on post-shock sound waves 
contributes substantially to this slowness. 

Attempts to extrapolate solutions with different viscous lengths must cope with the variation 
in the widths of transitions that represent shocks. These widths closely scale with the artificial- 
viscous lengths, but small deviations from exact scaling behavior may pose considerable barriers 
to extrapolation algorithms. These deviations from scaling are the main subject of Appendix C. 



2. Numerical Methods 



We refer to our two finite-difference schemes as "internal-energy" and "total-energy," after 
the underlying differential equations that they represent. We consider the one-dimensional Euler 
equations, written in the form 

dp dm 
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where m = pv is the momentum density (here per length); Cmt, the internal energy density; etot, 
the total energy density, and P, the pressure. The equation of state is given by the ideal gas law; 
generally we take 7 = 5/3. Equations (^) and @ are alternative expressions of the conservation 
of energy. In the internal-energy scheme, we apply equation (^); in the total-energy scheme, 
equation 

Since we want to capture shocks automatically and so avoid the computational difficulties of 
shock tracking, we use an explicit artificial viscosity to generate the entropy necessary for shocks. 
Following Norman and Winkler ( 1986| ), to the physical pressure P we add the term 



dv 



dv 

-Lic + Lr, mm ( — — , 
dx 



(5) 



where c is the adiabatic sound speed. The constants Li and L2 are viscous lengths for linear 
and quadratic artificial-viscous terms, and parametrize their strengths. Typically the artificial 
viscosity smears shock transitions over a region whose size is comparable to the viscous lengths. 
Conventionally the viscous lengths Li and L2 are chosen to be a small multiple of the gridsize Ax; 
we adopt this approach for our initial investigations. 

For a finite-difference representation of the equations of motion, we employ an operator-split 



scheme similar to the internal-energy scheme developed by Norman and Winkler ( Norman fc 



Winkler 1986| , ^tone &: Norman 19921) . They use a staggered spatial grid, defining scalar and 
vector quantities at different gridpoints; whereas our grid is unstaggered. 

The operator splitting proceeds as follows: for the parts of the equations (||-@) that represent 
advection, 

dip d[vtp] 

dt dx ' 



(6) 



we use the monotonic advection scheme of van Leer ( |1977| ). The rest of the terms are treated in 
three further substeps: pressure acceleration, artificial viscosity, and compressional heating. In 
each full timestep, we account first for each of these three terms (in the order specified), and then 
for the advection terms. Each substep generates new values of the physical variables, which then 
are used in the update for the next substep. 



Explicitly, then, the finite-difference equations take this form: first we represent the 
pressure-acceleration term by 



At 

2Ax 



(7) 
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To incorporate the artificial-viscosity terms we next use 

At 



■■-J --J 2A2; '-'^^^^ ~ ^^"^ 

{eint)j ^ {eint)j - -^^Qj ~ "i-i) (9) 

{etot)j ^ i(^tot)j - ^^{Qj+iVj+i - Qj-iVj-i), (10) 

where values of the momentum density — updated by pressure acceleration — enter into the 
right-hand sides through Vj = rrij/pj and Qj = pj{dv/dx)j[—LiCj -|- min (0, For 
the gradient term, we take {dv/dx)j = (fj+i — t'j_i)/2Ax. Then, in the compressional-heating 
substep, the energy density is changed by 

(emt)j 1-^ {(iint)j - -^^{eint)j{vj+l - Vj-l) (11) 

{etot)j ^ {etot)j - [{etot)j+iVj+i - {etot)j-iVj-i] , (12) 

where once again the latest values of {eint)j, {etot)j, and mj enter the right-hand sides. 

The resulting values of the physical quantities enter the advection substep. The van Leer 
procedure has the flux-conserving form 

V'j ^ v-j - x:^(-^i+i/2 - ^j-i/2)- (13) 

The fluxes are calculated in the upwind fashion 



Vj+1/2 
Vj+1/2 



V'i+i - ^ [vj+i/2^t/ Ax + Ij VV'j+i 



if < 0, 



(14) 



where ^^+1/2 = {vj+i +Vj)/2, and (as always) vj = rrij/pj. This advection scheme's monotonicity 
arises from van Leer's definition of the gradient-like term 

W - = I ^^^i+V2'^'^i-l/2/(V'i+l - i'j-l) if '^V'j+l/2#i-l/2 > 0, ^^^^ 

1 otherwise, 

where 5V'j+i/2 = ^j+i ~ i'j- The spatial accuracy of the van Leer procedure is second-order, except 
at extrema, where the accuracy must be reduced to first order to guarantee monotonicity. After 
the advection substep, the update through the timestep At is complete. 

In regions where the continuum solution is smooth, such difference schemes are consistent 
with the differential equations ( ^od 198^ ). When discontinuities are present, however, it is not 
immediately obvious that truncation errors should vanish in the Ax, At — > limit. If we regard 
truncation errors as products of these small quantities and derivatives of the solution, then clearly 
it is possible that the truncation errors may remain finite even in the small-gridsize, small-timestep 
limit. 
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Lax and Wendroff ( 19601 ) proved that single-step conservative difference schemes, if stable 



and convergent, converge to weak solutions of the corresponding systems of conservation laws. 
This theorem can be directly extended to include operator-split, multi-step conservative schemes 
as well (see Appendix A). We know of no such theorem or extension for non-conservative schemes, 
however. Below we present results that show that some non-conservative schemes (namely, our 
internal-energy type) do not converge to weak solutions of conservation laws. 

The operator-splitting method limits the overall time accuracy to first order, and the 
spatial accuracy of the substeps is (almost everywhere) second order. For smooth flows, then, 
asymptotically the error per step is 0((At)^, At(Ax)^). 

In general, the need for numerical stability requires constraints on the 



timestep At ( Richtmyer fc Morton 1967 ). We impose two constraints. First, the timestep 



must everywhere satisfy the Courant-Friedrichs-Lewy condition. 



At < r-^. (16) 

~ \v\+c ^ ^ 



Second, since the artificial viscosity takes the form of diffusion, we enforce a "diffusion-limiting" 
condition 

A. < , (1, 

Lic — 2L2 min {dv/dx, 0) 

To ensure these conditions, before making each update we choose the timestep to be 



At = min 



CAx D{Ax 



\2 



(18) 



|?;|+c Lic — min (c^w/Sx, 0) 

where C < I and D < I are safety factors. All of the calculations shown in this paper 
useC = D = 0.9. 

Despite these stability constraints, we have found that our total-energy scheme is unstable in 
the limit Ax — > when Li/Ax = constant. In our experience, any given unstable calculation may 
be stabilized by increasing the linear artificial viscous length Li. However, we discourage the use 
of the total-energy scheme for convergence studies. Imposing a diffusion-limiting condition (17) 



means that increasing the artificial viscosity forces smaller timesteps and greater computational 
expense. We use the scheme in the context of this paper to emphasize that our conclusions 
regarding the internal-energy scheme are based on that scheme's non-conservative nature. 

In the calculations that we show in this paper, we observe no indications of instability. 
Further, we have duplicated most of the total-energy calculations with the (likewise conservative) 
classic Lax- Wendroff scheme ( Lax fc Wendroff 196C| ) in two-step form ( Richtmyer fc Morton 19"67| ). 



(Some Lax- Wendroff calculations required reduced timesteps for stability, achieved by using a 
smaller safety factor D.) These verify that our results regarding convergence in the presence of 
linear artificial viscosity depend on the viscosity, not on the underlying numerical scheme. 
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3. Nonlinear Sound Wave 



Since the internal-energy equation (|3|) is not in the form of a conservation law, finite-difference 
schemes based on it do not conserve total energy exactly. Hence, one test of its results is to 
compare the calculated total energy with the total energy in the initial conditions. Under some 
conditions this scheme produces incorrect results. 

For the first test problem we consider a steepening wave. The initial conditions superpose a 
large-amplitude wave on a flat background: 

v{x) = 6{x) 

P{x) = 3/5 + <5(a;) (19) 
p{x) = l + 5[x), 



where 



5[x) = i '^o'' ' 1^ < (20) 

I 0, otherwise. 



For the amplitude we take 6q = 0.2; and for the width parameter, A = 100/13. In the 
small-amplitude limit of the inviscid equations of motion, this waveform is a purely right-moving 
wave. Because of non-linear effects, it steepens and forms a shock on its leading edge. (The grid is 
large enough that waves from the perturbation do not reach the edges; hence no energy flows in 
or out the computational volume.) The viscous parameters are Li/Ax = 1/2 and L^j ^x = 1. We 
refine the gridsize Ax and compare the results. 

Figure 1 illustrates the behavior of the total-energy error AEfoiai at various times as a 
function of gridsize Ax. At early times {t < 20) the error decreases proportionally to Ax, as 
anticipated for this first-order accurate scheme. At late times {t > 28), however, at small values 
of Ax the energy error tends toward non-zero values. For these cases, increasing resolution leads 
not to greater accuracy, but rather toward some incorrect solution. The internal-energy scheme 
cannot be performing correctly at late times. 

The shock forms at t = 23 it 0.5, estimated as follows. We repeat the calculations beginning 
with the same initial conditions but without artificial viscosity. To forestall the onset of numerical 
instabilities that develop in the absence of linear artificial viscosity, the timestep safety factor C 
is reduced to 0.4. We use the characteristics moving to the right from each gridpoint at the 
speed V + c to find the time of first crossing, which corresponds to shock formation. Close 
examination of the t = 24 curve in Figure 1 suggests that it may deviate from the well-behaved 
curves for t < 20. Once the shock forms, the non- vanishing energy error becomes increasingly 
apparent. This experiment strongly suggests that the internal-energy scheme fails when the flow 
contains a shock. 



Before discussing this failure further, we introduce another measure of errors in computed 
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solutions. Define the Cauchy error for a solution computed with gridsize Ax by 



eAxbP]= J tpAx{x,T) -i)^^i2{x,T) dx, (21) 

grid 

where ipAxix, T) is some physical variable (at position x and time T) as computed with gridsize Ax. 
When the numerical calculations converge to the solution of the continuum system, this error 
measure should behave as 0((Ax)"), where a gives the leading-order behavior of the error. For 
smooth flows, a is simply the order of accuracy of the scheme. When discontinuities are present, 
however, a depends also on the changes in the computed transition widths i, which should be 
on the order of the viscous lengths Li and L2. (We discuss this issue further in Appendix C.) 
If the jump across the transition is A^, then (on dimensional grounds alone) the shock should 
contribute (up to numerical constants) to the Cauchy error. When the viscous lengths Li 
and L2 are taken to be proportional to the gridsize Ax, the contribution to the Cauchy error should 
also behave as 0(Ax). (Note, however, that this is a norm-dependent statement. For example, 
in the commonly used C2 norm the error at shock transitions should behave as 0((Ax)^/^).) 
When computations are sufficiently well-resolved so that solutions converge as expected — that is, 
e[V'] oc (Ax)", where a has its expected asymptotic value — then they can be regarded as quite 
close to the Ax^O result ( [Finn &: Hawley 1989 ). Note, however, that the Cauchy error exploits 



no knowledge of the exact continuum solution. In practice one often assumes that convergence at 
the expected rate implies convergence to the correct solution. 

Figure 2 shows that the Cauchy error mcB-surcs for all three densities ttz, and. Cini ) decrease 
linearly with decreasing gridsize Ax. This shows that convergence of the Cauchy error alone 
does not imply correct behavior of the numerical scheme. (Since the Cauchy errors for all three 
quantities behave similarly, in the rest of this paper we illustrate Cauchy errors for the momentum 
density only.) 

Aside from the total-energy error, we have no indications that the internal-energy scheme 
behaves incorrectly. No numerical instabilities appear, and the performance at pre-shock times 
is exactly as anticipated. The error itself is relatively small — at i = 40 the error in total energy 
approaches approximately 0.4% of the initial energy in the wave — and has likely gone unnoticed 
by many users of internal-energy schemes. Since the scheme appears stable and convergent, this 
finite error clearly shows that non-conservative schemes need not converge to solutions of the 
conservation laws when discontinuities are present. 

Motivated by this hypothesis, we make the following modification: Instead of scaling the 
viscous lengths Li and L2 with the gridsize Ax, so that Li/Ax, L2/AX = constants, we treat the 
viscous lengths themselves as constants. This approach changes the differential equations being 
solved. Viscous effects remain finite even as the gridsize Ax vanishes, and so shocks are replaced 
by smooth transitions with finite width. With this change, the internal-energy scheme behaves 
properly at all times. We expect errors proportional to (Ax)^, since the scheme's first-order error 
is entirely temporal (i.e., 0(At), not 0(Ax)), and, by the diffusion-limiting condition (17), for 
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sufficiently small Ax the timestep At is proportional to (Ax)^. 

We repeat the steepening- wave problem, now choosing Li = 1/4 and L2 = 1/2. In Figure 3 
we see that, for sufficiently small gridsize Ax, the total-energy error AE total tends toward zero 
as 0(Ax)^ at all times. Figure 4 shows that the Cauchy error e[m] also decreases as 0(Ax)^ 
at all times. The magnitudes of these errors are much smaller than their counterparts (at the 
same gridsize Ax) from the previous (Li/Ax, L2I Ax = constants) calculations. This suggests 
that these finite-gridsize results are quite close to the solutions of the (finite-viscosity) continuum 
system. Equivalently, it suggests that the Cauchy errors measured in the previous calculations are 
dominated by the change in viscous lengths Li and L2, not by effects of the non-zero gridsize. 

Although the preceding discussion has focussed on the failure of a particular internal-energy 
scheme with particular computational parameters, the phenomenon is more general. We 
have varied the size of the timestep safety factors C and D and the strengths of the viscous 
parameters Li/Ax and L2/AX (including setting L2 = 0) with no improvement in results. Further, 
we have implemented versions with Norman and Winkler's staggered grid, as well as several other 
modifications to the difference equations (Norman & Winkler 198(:, Stone &i Norman 1992| ), also 
without improvements in results. (Appendix B details the modifications we have tested.) 

It is possible that the scheme may converge toward the inviscid solution if we choose to vary 
the gridsize Ax and viscous lengths Li and L2 such that Ax,Li,L2^0 but Li/Ax, L2/Ax^oo. 
In this situation the width of the shock representation should decrease with increasing resolution, 
but the number of gridpoints in the representation should increase. We do not investigate the 
behavior of the scheme using this procedure. The slowly decreasing values of Li and L2 impose 
penalties on computational efficiency [via the diffusion-limiting timestep condition (17)], and it 
is unclear whether convergence toward a physically reasonable solution can be readily observed 
without choosing the viscous lengths to be nearly constants. 

To summarize, when the flow contains shocks (and the limiting equations are inviscid), the 
internal-energy scheme converges toward a solution that is demonstrably wrong. This failure 
arises because the solution contains discontinuities and the difference scheme is non-conservative. 
We can recover acceptable behavior by changing the equations so that the limiting solution does 
not contain discontinuities. In the next section we study another test problem in which the 
internal-energy scheme demonstrably fails, and also show that the total-energy scheme succeeds. 



4. Shock Tube 



Because we associate the failures of the internal-energy scheme with the presence of shocks, for 
the second test problem we consider a shock tube, as introduced by Sod ( 1978|) . At t < a static, 
impermeable interface separates two states of a gas at rest. At i = the barrier disappears, and a 
shock then forms. Because the problem admits a similarity solution that can be found exactly, it 
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is an ideal test problem for shocked flows. We test the performance of both the non-conservative, 
internal-energy scheme and the conservative, total-energy scheme. 



The exact solution can be constructed as follows (Landau & Lifshitz 1987): Since neither 
the equations of motion nor the initial conditions contain a characteristic length, all interfaces 
must travel at constant speeds away from the initial discontinuity. A shock propagates into the 
low-pressure region; a rarefaction wave, into the high-pressure region. The gas between these 
two transitions travels at a constant speed and is uniform, except at a single point at which the 
density may be discontinuous, the contact discontinuity. The Rankine-Hugoniot jump conditions 
relate the upstream velocity, pressure, and density to the uniform middle region's conditions. The 
transition across the rarefaction wave can be obtained by noting that the flow in this region is 
adiabatic, and so the physical quantities can be determined with the method of characteristics. 
Only one possible shock and one possible rarefaction wave allow the velocity and pressure to 
match in the middle region; we determine these numerically, since the relations are implicit. 

With the exact solution in hand, we can test both internal-energy and total-energy schemes. 
We use Sod's initial conditions. On the left side of the grid, the pressure and density are given the 
values 1; on the right, the pressure is 1/10 and the density is 1/8. (At the interface we give these 
quantities their averaged values.) In this section only the equation of state is given by a 7 = 7/5 
ideal gas. 

Figure 5 dramatically illustrates the failure of the internal-energy scheme. (In the shock-tube 
calculations, we use Li/Ax = L2/AX = 3/2.) We show the density in the region around the shock 
jump at t = 1 and compare the results of the two numerical schemes with the exact solution. 
Although the computed solutions appear well-behaved and free of instabilities, two problems with 
the internal-energy solution are apparent: first, the shock transition is in the wrong place (the 
shock velocity is 0.2% too low), and second, the post-shock density is wrong (0.3% too high). 
Further, the total-energy scheme gets both of these right; aside from a small post-shock oscillation, 
it exhibits no numerical difficulties. In Figure 6 we show the total-energy error accrued by the 
internal-energy scheme at various resolutions, and again this diagnostic verifies that the scheme 
does not converge toward a correct solution. 

In Figure 7 we show the Cauchy error measures eAa;["^] associated with both schemes' 
calculations. For each, the errors decrease with decreasing gridsize; in fact, the quantities are 
almost identical. As before, however, we already know that the internal-energy scheme does not 
converge to the correct solution. 

The rate of convergence appears poorer than 0(Ax). We trace this slow convergence to the 
presence of contact discontinuities, which can present special problems to numerical schemes. 
We observe the correct rate of convergence when we alter the initial conditions to remove the 
discontinuity. This change removes much of the shock-tube problem's usefulness as a test problem, 
though, since the time evolution is no longer given by a similarity solution. (In the next section, 
in which we discuss rates of convergence, we use a test problem with smooth initial conditions 
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that does not contain contact discontinuities.) 

Unlike shocks, contact discontinuities do not develop in smooth flow. Similarly, the smoothed 
regions that represent contact discontinuities do not tend to steepen. Thus effects of non-zero 
gridsize cause them to spread without limit, rather than toward some steady state. For a contact 
discontinuity in an isobaric medium, translating at a constant speed, no length or time scales enter 
except the width w of the representation and the gridsize Ax. Thus a scaling holds among w, Ax, 
and the time t: we may write that w{Ax,2t) = 2w{Ax/2,t). It is clear that if the representation 
spreads without limit, then at all times w{Ax,2t) > w{Ax,t) and so w{Ax/2,t) > w{Ax,t)/2. 
Since the Cauchy error accrued at these discontinuities is proportional to their widths, the 
Cauchy error cannot scale linearly with the gridsize Ax. (We might expect that this diffusion 
process would lead the width w to behave asymptotically [as t^oo] as t^^'^. Instead, however, 
using the total-energy scheme we find that a calculation involving only the advection of contact 
discontinuities suggests that asymptotically w varies as t^^^.) 

For the shock tube, as in the steepening wave, the internal-energy scheme appears stable 
and convergent but fails to converge toward a correct solution. The Lax-Wendroff theorem shows 
that the conservative, total-energy scheme should converge toward a correct solution, and all 
indications are that it does so. 

In brief, when the artificial-viscous lengths Li and L2 are proportional to the gridsize Ax, 
the calculations of our internal-energy scheme do not converge to the known shock-tube solution; 
while those of our total-energy scheme do approach the exact solution. This result reinforces 
our finding that the internal-energy scheme fails in the presence of shocks. Further, the Cauchy 
errors (^l]) for the (correct) total-energy scheme and the (errant) internal-energy scheme are nearly 
indistinguishable. Use of convergence tests alone cannot guarantee convergence to the correct 
solution. 



5. Effects of Diffusion on Convergence 

Since the internal-energy scheme is unreliable when shocks are present, we employ the 
total-energy scheme exclusively for our discussion of convergence rates. The diffusive effects of 
linear artificial viscosity introduces a noticeable amount of diffusion even to linear perturbations. 
We first discuss in some detail the impact on small, linear waves. Then we show the effect on 
convergence rates for a wall shock — a test problem containing both a shock and small-amplitude 
waves. 



5.1. Linear examples 
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The presence of linear artificial viscosity can have a surprisingly large effect on the observed 
convergence rate. To illustrate these effects, we consider small-amplitude waves. First we exhibit 
convergence rates for calculations of such waves, and then show consistency between rates in this 
calculation and rates derived using an analytic treatment. Finally we comment on the viscous 
lengths necessary to yield convergence rates close to the expected asymptotic result. 



For the numerical calculations, we take as initial conditions the right-moving wave (1!; 
discussed above in section 3, with amplitude 6o = 0.05 and width parameter A = 0.4. We calculate 
to time t = 0.6 using the total-energy scheme with linear artificial viscosity only — the quadratic 
artificial viscous length L2 is 0. Figure 8 shows the Cauchy errors e[m]. Note that instead of 
plotting eAxlw-] as a function of Ax (as in Figures 2, 4, and 7), here we show it as a function 
of the viscous length Li. The computations use Li/Ax = 4, chosen so that the resolution is 
sufficiently high that purely numerical errors (as compared with infinitely resolved computations 
with those values of the viscous length Li) are negligible. Note that as Li decreases the decrease 
of the errors tends toward the expected linear rate (cf. the line segment corresponding to 
0{Ax) dependence), but only rather slowly. To quantify this tendency, for triplets of calculations 
we define (3{Li) = log(e2Li ["ij/e^j [m])/ log 2, which approximates the logarithmic derivative of e^^ 
with respect to Li. In the asymptotic regime, (3 approaches the order of accuracy. We show /3 as a 
function of Li in Figure 9. Note that even at the smallest value of Li, {3 is still rather far from the 
expected value of 1 even though almost two orders of magnitude separate Li from the problem's 
natural length scale A. 

A simple analytic model shows such convergence rates are a consequence of the diffusive 
effects of linear artificial viscosity. Consider the equations of motion (||-|^ linearized about a 
constant background, which (in terms of the small quantities pi, vi, and Pi) read 

dpi ^ _y^^_p^^ (22) 

dt dx dx 

dvi _ dvi I dPi ^ ^ ^ d'^vi 

dt ^ dx po dx ^ ^ dx"^ 

dPi ^dvi dP, 

(The background is specified by the constants vq, po, Pq, and cq = (7^0/^0)"^'^^- Note that Li is 
the linear viscous length, not a perturbation quantity.) We include the artificial-viscous terms, 
but to leading order in the small perturbations the quadratic artificial viscosity vanishes. In the 
absence of linear artificial viscosity (Li = 0) the linearized equations take the familiar Euler form, 
and the three coupled equations can be rewritten as three decoupled equations, 

-df = 5^ 

9X2,3 I , nC^X2,3 

= 

where Xi = Pi ~ -Pi/cq and X2,3 = {Pi/cq ± pqVi)/{2co), for modes that propagate along the 
three characteristic paths, at speeds vq and vq it cq. In the presence of linear artificial viscosity. 
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performing the analogous transformation leads to three still coupled equations: 

-m = 9^ 

^ = -(^o + ^o)9^ + ^5^(^2-X3) (28) 

^ = -^'''-''^^ + ^d^2ixs-X2). (29) 

Note that xi) the mode that remains still with respect to the background velocity vq, is unaffected 
by artificial viscosity. We next assume that (as in the wave example above) one of the remaining 
modes is sufficiently larger than the other two so that we can neglect their influence on its 
propagation. Calling it simply we thus obtain the equation of motion 

^ + (.o±co)^ = ^^. (30) 

The left-hand side of this equation describes just the advection of the quantity X) so that (with 
the identification of the diffusion constant v = LiCq/2) the equation is simply the heat diffusion 
equation 

for the quantity xiUii) i^i the frame of reference that moves at the speed vq it cq. 

We study the convergence of solutions of the heat equation (31) as the diffusion coefficient v 



tends to zero. We consider the time evolution of a Gaussian pulse. An exact solution for t > is 

X{y,t;u) = A (^^;^)'^' e-yy^^'^'^^ll (32) 

where A is the initial (t = 0) amplitude and Aq is the initial width of the pulse. With this solution 
it is straightforward to calculate the Cauchy error (at a fixed time t) associated with a sequence 
of solutions in which the viscous parameter varies. In Figure 10 we show the Cauchy error €u[x] 
evaluated at t = 0.6 using an amplitude A = 0.05 and width Aq = 0.3. (We choose different values 
for the width parameters in this treatment and in the calculations described above because the 
shape of the initial conditions differs.) Again the convergence rate /3, shown in Figure 11, becomes 
close to the expected value of 1 only at quite small values of the diffusion u. In the context of the 
flat background of the previous calculations, we require values of the viscous length Li = Iv jcQ 
much smaller than the natural length scale Aq to achieve values of /? close to 1. 

To determine how small the diffusion coefficient v — or, equivalently, the viscous length L\ — 
must be in order to achieve a convergence rate /? close to the asymptotic value of 1, we examine 
the pointwise convergence of the exact solution ( ]3^ ) of the linearized problem. That is, for any 
values of {y,t), consider the behavior of xiu^'^j^) as i/— >0. We assume that, for sufficiently small 
values of u, we may approximate 

ov 2 ov 



X{y,t;u) - xiy,t;0) - i^^{y,t;0) + —^iy,t;0). (33) 
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We determine the value oi u = v{y, t) at which the first and second terms on the right-hand 
side are equal. The asymptotic regime for convergence lies at values of much less than V. The 
calculation itself is elementary and yields 



^0 



2y' 



^0 



V-12y2A2 + 3A4 



(34) 



The analogous quantity for the viscous length is 



Li(2/,t) = Ao^ 



Cot 



2„,2 



2A4 



V - 12y2A2 + 3A4 



(35) 



Loosely, the convergence rate for the Cauchy error is a weighted (spatial) average of the 
pointwise convergence rate. Since the main contribution to the error comes from a region within 
a few widths Aq of the spatial origin, we replace the dimensionless ratio in equation ( |35| ) by 1/2. 
Thus, the asymptotic regime for the Cauchy error should lie at values of the viscous length Li 
much less than 

17 (36) 
2 cot ^ ^ 

Note that smaller values of Li are required as the time t increases to achieve the same convergence 

rate. 

In practice, one may wish to increase the resolution of a computation until the convergence 
rate P is within some tolerance e of the asymptotic value (here, 1). Assume [as above, in 
equation (^3|)] that for small L the error behaves as 



eL (xLL±L^ 



(37) 



where L specifies the value of L at which the two terms on the right-hand side contribute equally. 
The logarithmic derivative of the error ei with respect to L is then simply 



L dcL _ 1±2L/L 
e^~dL ~ 1±L/L 



1±L/L. 



(38) 



To attain a convergence rate within e of 1, we must use values of L less than approximately eL. 
(Comparison with the results shown in Figures 9 and 11 suggest that this estimate is about right 
for the analytic model and too high for the numerical calculations by a factor of a few.) 



5.2. Wall Shock 



As noted above, even for the reliable total-energy scheme our results for the shock-tube 
problem do not yield the expected 0{Ax) rate of convergence. This is a consequence of the shock 
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tube's discontinuous initial conditions. For our final test problem, we choose one that contains a 
shock, smooth initial conditions, and small-amplitude waves: a wall shock. 

Gas impinges on a wall at x = 0. For uniform, supersonic flow from the right (x > 0), a 
shock will form and move outward. For initial conditions we place a transition between density, 
momentum, and energy values appropriate to a Mach 2 steady-state shock. (In the equation of 
state we again use 7 = 5/3, as in the study of the nonlinear sound wave.) We mediate between 
upstream and downstream states with a smooth transition function: 

v{x) = /(0,-9/8;x) (39) 

P{x) = /(57/20,3/5;x) (40) 

p{x) = /(16/7,l;x), (41) 

where the transition function is 

{a, if X < 3 

(a + 6)/2 + [(6-a)/2]tanhf|(x-6), if 3 < x < 9 (42) 
b, if X > 9. 

At the left edge (at x = 0) we use reflecting boundary conditions, but for the tests reported here 
we stop the calculations at T = 2 before any physical perturbations reach either edge of the 
computational grid. 

Although the initial conditions do not represent the exact, steady-state profile for a shock 
transition, one quickly forms and translates away from the wall. Waves form and propagate away 
from the shock toward the wall. Figure 12 shows these structures in the density. 

We have performed many calculations parameterized by the viscous length Li = L2 = L 
and the gridsize Ax. We use these to form two types of sequences: first, with the viscous 
parameter L/Ax held constant, and second, with the viscous length L held constant as Ax varies. 

Figure 13 shows the Cauchy error eAx•[^T^], again as a function of the viscous length L. Each 
of the five solid lines represents a sequence with constant L/Ax (= 1,2,4,8, 16). First, note that 
over the range of gridsizes Ax in which we can reasonably calculate, the convergence rate tends 
toward 0(Ax) but does not approach it unambiguously. Note also that the magnitude of the 
errors is almost entirely a function of the viscous length L — the solid lines overlap. That is, the 
particular choice of gridsize has only a small effect on the Cauchy error. The difficulty of achieving 
the expected convergence rate in these sequences must be an effect of the changing values of the 
viscous lengths, not of the finite resolution. 

We attempt to isolate the effect of the shock transition on the Cauchy error by changing the 
limits of integration in its definition (^) to include only the transition region or to exclude it. 
Figure 13 also shows in dotted lines the error associated with the transition region (x > 7), and 
in dashed lines with the transients in the rest of the flow (x < 7). Note that the transition region 
yields the expected 0(Ax) convergence rate at relatively large values of the gridsize Ax. The 
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more slowly decreasing contribution to the Cauchy error appears in the region of the transients. 
The Cauchy error for the shock transition should asymptotically be proportional to the viscous 
length L, as discussed above, and so its convergence rate implies that transition widths are indeed 
scaling with L. (Appendix C discusses the small deviations from exact scaling in shock transition 
regions.) 

Figure 14 shows again Cauchy errors as a function of gridsize lS.x. Here, however, each solid 
line represents a set of calculations with a fixed value of the viscous length L. We show only 
calculations with the gridsize Ax sufficiently small that the timestep At is always determined by 
the diffusion condition ([T7|) ; then the errors should (asymptotically) be 0((Aa;)^). It is clear that 
these calculations lie in the asymptotic regime: the errors are nearly proportional to (Ax)^. 

We use the same calculations to construct Figures 13 and 14. When considered as part 
of sequences with constant L, they unambiguously lie in the asymptotic regime; whereas when 
considered as part of sequences with constant L/Ax they do not. We illustrate this dichotomy 
in Figure 15: each circle represents a calculation, and lines represent the sequences to which 
each belongs. For triplets of calculations in sequences, we calculate the numerical convergence 
rate /3(Ax) = log(e2Aa:["i]/eAx["i])/log2, which approaches the order of accuracy in the 
asymptotic regime. Values of /3(Ax) are shown on the lines that define the appropriate sequences 
of calculations. For each constant-L sequence, where we expect the errors to behave asymptotically 
as 0((Ax)^), we find that /? = 2 to several decimal places. For sequences of L/Ax = constant 
calculations the closest we can achieve to the expected value (/3 = 1) is /3 = 0.93. 

The Cauchy errors in the L/Ax = constant sequences are larger in magnitude than those in 
the constant-L sequences (compare Figures 13 and 14). The smallness of the latter errors and 
the attainment of the expected convergence rate shows that errors in the former sequences are 
dominated by the effects of changing viscous lengths, not of insufficient resolution. 

The apparently poor convergence rate in the post-shock region can be related (at least 
in part) to the discussion of small-amplitude waves above in section 5.1. The features in the 
post-shock region (cf. Figure 12) consist of relatively small-amplitude waves. A considerable 
contribution to the integrand defining the Cauchy error (^) comes from the region near x = 3, 
where a left-moving wave lies. From the previous discussion, for this wave alone we require much 
smaller viscous lengths L to attain a numerical convergence rate (3 within a few percent of 1. (We 
have deemed this to be computationally infeasible.) 

This explanation implicitly assumes that the waves generated as the shock forms do not vary 
with the viscous length L. The width of the shock transition, varies as the viscous length, however, 
and so some post-shock waves (also varying as 0{L)) must arise from the different transition 
widths. Along a sequence in which L/Ax = constant, these waves with short length scale should 
make a contribution of order 0(Ax) to the Cauchy error. The effect of the other (long length 
scale) waves on the numerical convergence rate continues to increase with time (i.e., fi is further 
from 1), and so ultimately the effect of these should be greater than that of the short- wavelength 
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waves. 

By studying this relatively simple test problem, we find that we can readily attain the 
expected 0((Aa;)^) convergence rate in sequences in which the viscous lengths Li and L2 are 
constants. We emphasize that in such sequences the Ax— >0 limit does not correspond to the Euler 
equations In sequences in which the ratios Li/Ax and L2/AX are constants, attaining the 

expected 0{Ax) convergence rate appears more difficult. The difficulty is tied to the changes in 
viscous lengths, not to numerical errors introduced by non-zero gridsizes. We attribute it to the 
diffusive effects of linear artificial viscosity on waves that are generated as the shock forms. 



6. Conclusions 

We have investigated the use of convergence tests to assess the accuracy of finite-difference 
solutions to the Euler equations in the presence of shocks. We find that a non-conservative 
internal-energy scheme does not converge to solutions of the continuum equations. Use of the 
Cauchy error does not detect this significant failure. The Lax-Wendroff theorem applies to our 
total-energy scheme (as shown in Appendix A), however, and so it provides reliable results. 

Our identification of this failure with the presence of shocks leads us to introduce a 
modification of the equations of motion, in which the artificial-viscous terms are held constant 
as the gridsize varies, rather than scaled with the gridsize. In such modified systems, smooth 
transitions replace shock discontinuities. With the discontinuities removed, the internal-energy 
scheme converges as expected, without the indications of the problems that signal the failure of 
the unmodified system. 

For well-behaved schemes, convergence at the expected rate can be used as a clear indicator 
of accuracy, a particularly unambiguous criterion in automated contexts. Diffusion due to linear 
artificial viscosity can lead to convergence difficulties. By investigating a wall shock with smooth 
initial conditions, however, using our definition of Cauchy error to provide the criterion for 
convergence, we readily achieve the expected convergence rate when the artificial-viscous terms 
are held constant, and approach it when the artificial-viscous terms vanish with the gridsize. 

Finally, we have investigated the possibility of using the calculations that establish our 
numerical convergence rate in extrapolations. In sequences in which the artificial viscosity scales 
with the gridsize, we find that transitions that represent shocks deviate in several ways from exact 
scaling. Our findings on extrapolations and scaling can be found in Appendix C. 
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A. Lax-Wendroff Theorem, Extended 



We discuss in this appendix the extension of the Lax-Wendroff theorem to muhi-step, 
operator-spht schemes. This theorem was originally presented by Lax and Wendroff ( 19601 ); our 
discussion follows their treatment. 

Write the system of conservation laws as 

^ = ff' ^^^^ 

where u is an unknown vector function of x and t and / is a given vector function of u. The 
theorem requires that the difference equations be written in the form 

v(x,t + At) = v(x,t) + -^Ag, (A2) 

Ax 

where A;; = g{x + Ax/2) — g{x — Ax/2). The numerical flux g is determined in the following 
fashion: g{x + Ax/2) = f_/+2) • • • ,vi), where G is a vector-valued function of 21 vector 

arguments, which is related to / by the single requirement that G approach /, in the sense that 

lim G{vi,V2,...,V2i) = f{v). (A3) 

all Vj^v 

If, in the limit Ax, At — > 0, v{x,t) converges (boundedly, almost everywhere) to some 
function u{x,t), then u{x,t) is a weak solution of (|A1|). 



We claim that the proof also holds if the numerical flux g is given by a vector-valued 
function G of vector arguments and a (scalar) parameter, A = At/ Ax, such that 

lim G{vi,V2,... ,V2r,X) = f{v). (A4) 

^t),A— +const. 

Operator-split schemes satisfy this condition, as we now explain. Write the continuum 



equations (Al) as 



du _ df^ 
dt dx dx 

where we take the number of substeps to be 2 (for notational simplicity). The difference equations 
will consist of two substeps to represent the equations 

dt ~ dx ^^^> 

sequentially. That is, the first substep takes the form 

= V, + ^^Ag^^ {{v,}) , (A8) 
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where {vj} represents the set of all relevant vj. The resulting set {vj} enters the second substep 

(A9) 



We assume that the substeps' numerical fluxes g^^^ and g^'^^ are consistent with the substep 
fluxes f^^^ and /^^^ in the sense of (A3). 



The split scheme can be rewritten in the awkward one-step form 



Vi + —Ag^^^ {{vj}) + —Ag 



Ax 

Vj + —A 
^ Ax 

A^ 

Vj + -;—Ag, 



' Ax 
5« ({.,}) + <7(^) 



(2) 



V, + ^Ag(^H{v,}) 



+ ^A,(i) ({.,}) 



^ ' Ax 

where g{{vj}, At/ Ax) is defined by the quantity in the square brackets. 

Since the scheme can be written in the conservative form (|A2| ), we now establish that 



(AlO) 
(All) 
(A12) 



the numerical flux Q satisfies the extended consistency requirement (A4). Let all of its vector 
arguments approach the value v, while At/ Ax approaches some constant, finite value. Then Ag^^^ 
must vanish (since g^^'^ is a continuous function of its vector arguments where all of them are 
equal), and so essentially G becomes g^^^ {{vj}) + 5*-^^ {{vj})- Since each of g^^'^ and g^'^^ are 
consistent with f^^^ and f^'^\ then G must be consistent with /. 



B. Internal-Energy Scheme, Variations 



In this appendix we describe changes made to the internal-energy scheme, whose basic form is 
given in section 2. We emphasize that none of these changes affects the conclusions of sections 3 
and 4: the scheme still converges to solutions that are demonstrably incorrect. 

1. The viscous lengths Li and L2 and the timestep safety factors C and D were varied in ranges 
given. 

parameter range 
~T[ 0.1-1 
L2 1-32 
C 0.3-0.9 
D 0.008-0.9 

2. A staggered grid allows cleaner implementation of the van Leer advection scheme. The vector 
quantities — velocity and momentum — are stored at gridpoints midway between the gridpoints 
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for the remaining (scalar) quantities. This requires that every difference equation be modified: 
Pressure acceleration (|^ becomes 



At 



1714 



The artificial-viscosity difference equations dl, 0) become 



At 



m 



,+1/2 ^ mj+i/2- —{Qj+i-Qj 



{^int ) j 



Ax 
At 



(B13) 

(B14) 
(B15) 



where now {dv/dx)j = (fj+1/2 



12)1 ^x. The compressional heating equation (11) becomes 



At 



(B16) 



The van Leer advection procedure is virtually unchanged for scalar quantities. The velocity at the 
staggered gridpoints is taken to be 'fj+1/2 = '^f^j+i/2l {Pj + Pj+i)- For the momentum (a vector 
quantity), the velocity at the other gridpoints is taken to be Vj = (mj_i/2 + "^j+i/2)/(2pj)- 



3. As in Norman and Winkler ( 19861) and Stone and Norman ( 1992 ), the pressure-acceleration 
term in the momentum equation (0) can be replaced by 



dv 
dt 



1 dP 

p dx 



The difference equation for the physical pressure ( pi3 ) must be replaced by 
and that for the artificial- viscosity term ( [B14| ), by 



2At P,+i-P, , 

Ax pj+i + pj 



2At Qo+i — Qj 



(B17) 



(B18) 



(B19) 



4. In Norman and Winkler ( 1986| ) and Stone and Norman ( 1992| ), energy conservation is improved 
by using a "time-centered pressure" (as in P"+i/2 = (p" _|_ P"-+i)/2) on the right-hand side of 
the compressional- heating term of the internal-energy equation (^), then applying the equation of 
state to replace the pressure with the internal-energy density. The substep then can be expressed 

as 

, . l-At(7-l)(z;,+i/2-t^,-i/2)/2Ax 

which replaces (pH). 



(B20) 
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5. Consistent-transport advection ( Norman, Wilson fc Barton 1980| ) is said to improve local 
conservation of certain quantities. The advection of mass density remains the same. The remaining 
quantities are treated "consistently" with the mass density. In the van Leer advection substep, 
the momentum, rather than the velocity, moves the quantities over the grid. In the procedure 
for calculating the fluxes the momentum density field replaces the velocity field, and (to 

compensate for the extra factor of mass density) the specific densities ^//? replace the densities if) 
themselves. 



C. Extrapolation and Transition-Region Scaling 



In the main body of this paper, we adopt the attitude that we may regard a calculation as 
reliable if we can show that it lies within the asymptotic regime (i.e., if its errors are converging 
as expected). Since such a determination uses several calculations at different gridsizes, it is then 
straightforward to attempt to calculate an even more accurate solution by extrapolating the set of 
calculations. 

We first consider the wall-shock problem of section 5.2. Because its physical structure is 
relatively simple — one shock wave and detached waves — we investigate the possibility of using the 
many calculations to perform extrapolations to the Ax— >0 limit. Where the flow is smooth, an 
extrapolation in powers of the gridsize lS.x is valid. If the calculations extend into the asymptotic 
regime, then the extrapolations should be quite accurate. 



C.l. Extrapolation at constant L 



First, consider extrapolations of sequences in which the viscous lengths Li and L2 are 
held constant. We use calculations for which timesteps At are proportional to (Ax)^; then the 
truncation error associated with each is of order 0((Ax)^). We assume that, at each gridpoint in 
the best-resolved calculation, the results can be expressed as 



i^^x{x, T) = ^{x, T) + a2(Ax)" + a3(Ax)3 + • • • , (C21) 



with some coefficients (ifi — (x,T) ( [Richardson fc Gaunt 19271 ). (The linear term is absent 



because the scheme is second-order.) Extrapolation of N calculations finds the value of 
^{x,T) = ipAx=o{x,T) by eliminating the first — 1 error terms. Where values are required 
from calculations at points other than gridpoints, we interpolate using polynomials to an order of 
accuracy at least that of the eventual extrapolated result. The extrapolations are well-behaved 
and apparently more accurate than their component calculations. We illustrate this accuracy in 
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Figure 16 by means of a Cauchy error obtained using extrapolations. The extrapolations, for 
which Li = L2 = L = 1/8, use calculations (with gridsizes differing by factors of 2) from Ax = 1/16 
down to the value of the gridsize indicated on the horizontal axis. The error decreases markedly 
with increasing resolution (which also corresponds to more extensive extrapolation). 



C.2. Extrapolation at constant L/Ax 

Next, consider extrapolations in which the viscous parameter Li/Ax = L2/AX = L/Ax is 
held constant. (For these extrapolations we assume that the errors take a form similar to that in 
equation (|C21| ), but include the term linear in the gridsize Ax, because of the first-order accuracy.) 
Figure 17 shows the density as extrapolated from a sequence with L/Ax = 2. (This sequence 
includes the calculation illustrated above in Figure 12.) Qualitatively, the results appear poor in 
the region of the shock, but sensible in the remainder of the flow. The sharp features of the waves 
are ungarbled by the extrapolation, even though some features must result from the formation of 
the shock discontinuity. In the shock region, representations of shocks approach discontinuities in 
the Ax— >0 limit. Hence, although the Cauchy error measure eAa;[^] for the entire flow approaches 
the expected 0(Ax) behavior, near the shock the pointwise behavior of the error cannot lie in the 
asymptotic regime. 



C.3. Transitions in shock regions 

Although pointwise extrapolation of a sequence of calculations fails near shocks when the 
viscous lengths Li and L2 change within the sequence, one might try to devise a special procedure 
for extrapolation in shock regions. Any such method must consider the expected scaling of the 
transition region with the viscous lengths. We investigate the mechanisms that alter this scaling. 
We begin with a steady-state problem, the shock tube, to investigate short length scales, then 
move to the wall-shock problem for the influence of large-scale gradients. 

Consider a steady-state shock with uniform upstream and downstream conditions. The 
shock travels at constant speed, and the transition width for the shock representation must be 
determined by the viscous lengths Li and -L2. For simplicity, we take Li = L2 = L. Transitions 
calculated at different viscous lengths L then are related by a scaling relation: 

V^(x,i;L) = ^(^^|^) (C22) 
in the region of the shock, located at Xs, where the function ^ is independent of L. 



- 23 - 



Given two calculations \(j^^\x,t) and ^l^^'^\x,t) made with different gridsizes Ax, viscous 
lengths L, or both, we measure the deviation from exact scaling by means of the function 

1^(1.2) [V,] (e, t) = (x« + L(i)e, t) - ^P^ll,, {xf^ + L('k, t) (C23) 

in the vicinity of ^ = 0. Calculating this deviation function requires knowing the shock 
positions xi^^ and xi^^ In practice, when comparing computations with different viscous lengths L, 
we assume that xi^^ = xi^^ and approximate the shock position by locating the intersection of the 
two shock representations. Determining Xg in other ways does not change our results significantly. 



The deviation function vanishes if the calculations obey the exact scaling relation ( C22 ), since 
then V(i)(xi^^ + LW^,i) = iP^^\xP + L^'^k,t) =V'(0- When it is more convenient to work with a 
scalar deviation measure, we calculate a norm of 

D(1.2)[^] = r ^(1,2) ^(^24) 



To limit the investigation to the shock region, the parameter r] should not be chosen too large; we 
take 7] = 5. 

We begin our discussion of deviations from scaling with the shock-tube problem, in which 
the shock is truly steady-state. Our comparisons of calculations fall into two categories: first, 
we compare calculations with the same viscous length L; and, second, we compare calculations 
with different values of L. In the second type, comparisons between calculations with the same 
value of the ratio L/Ax are a useful special case. (For illustrations we focus on the sequence 
with L/Ax = 4.) 

Before describing the causes of scaling deviations, we first discuss some useful properties 
of shock-tube calculations. In this problem the only length scale set by initial conditions and 
equations of motion is the artificial- viscous length L. Thus solutions of the continuum equations 
(including the effects of artificial viscosity) obey the scaling property 

V'(x, t; L) = ijj{ax, at; aL), (C25) 

with a arbitrary, where the spatial origin is chosen at the point of the initial discontinuity. 
The corresponding finite-difference solutions have an analogous property. A difference 
solution ijj'j = tpAxU'^Xjtn; L), which approximates the continuum solution ipiijAxjtn), may 
equally well be interpreted as the result of a calculation in which the values of the viscous length L, 
the gridsize Ax, and the timesteps At are all scaled by a: ipj" = ipaAxiJcf'Axjatn'jCrL). Hence, 
two calculations with the same value of the viscous parameter L/Ax made to the same time T 
at different gridsizes (Ax)i and (Ax)2 are completely equivalent to two calculations with the 
same Lj Ax and (Ax)i, but made to the times T and ((Ax)i/(Ax)2)T, respectively. 



Truncation error (shock tube) 
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We have identified several causes of deviations from exact scaling. First, if the two viscous 
lengths entering into D[ip] are the same, then actually it does not measure deviations from 
scaling. Instead, it measures the truncation errors proportional to (Ax)^ of computed solutions 
(as compared to the "exact" solution with that same finite value of L). 

We estimate the magnitude of these truncation errors in a calculation at a given gridsize Ax 
by differencing it with a calculation at the smaller gridsize Ax/2 (with the same value of the 
viscous length L). In Figure 18 we show this estimate of truncation errors as a function of the 
scaled spatial coordinate ^ = (x — in the shock-tube problem for calculations (to the fixed 

time T = 1) at gridsizes Ax = 1/100, 1/200, 1/400, 1/800, all with the viscous length L chosen 
so that L/Ax = 4. Note that the results are all quite similar despite the different values of L 
and Ax. Along a sequence in which Lj Ax = constant, the effect of truncation error is a constant 
in the scaled coordinate ^. 

We explain this result as follows: Interpret the results as a sequence in time of truncation 
errors at a single gridsize and viscous length. If this chosen gridsize is Ax = 1/800, then Figure 18 
shows results at times t = 1/8, 1/4, 1/2, and 1, respectively. It then indicates that when the 
transition is near steady-state, the truncation error remains nearly constant in time. 

Time- dependent relaxation (shock tube) 

The solid line in Figure 19 shows the deviation function D[p\ that compares transitions 
calculated with viscous lengths L^^^ = 1/25 and L^^-* = 1/50, with the gridsize chosen sufficiently 
small (Ax = 1/1600) so that the 0((Ax)^) truncation errors are small. This illustrates a second 
source of deviations from exact scaling. We identify it as an approach to steady state, as the 
discontinuous initial conditions relax toward a steady-state computed transition. (The shape of 
the functions in Figure 19 show that transitions are steeper than the idealized, scaling case.) This 
relaxation is a consequence of the artificial viscosity, and so it is reasonable that its time scale 
must be on the order of L/u (where u is some characteristic velocity), which is the only time 
scale available in the problem. (Alternatively, we may argue that the diffusion time scale for a 
length i is on the order oi £'^/T>, with diffusion coefficient V oc uL. Since the computed transition 
width evolves from zero to a value on the order of L, the resulting time scale is on the order 
of L/u.) When the ratio L/T is smaller, the transition should be closer to steady state [i.e., the 
exact-scaling transition function (1C22| )]. 

Figure 19 also shows (dashed line) the deviation function calculated with the same values of 
the viscous length L but with larger gridsizes, chosen with the same value of the ratio L/ Ax = 4. 
With such large gridsizes, the 0((Ax)^) truncation error associated with either calculation is 
much larger in magnitude than the deviation function. The deviation function, however, is nearly 
identical to that in the more accurate, small-gridsize calculations. This shows a consequence of the 
similarity in truncation errors as a function of ^ = (x — Xs)/L along sequences with fixed Lj Ax: 
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Write the result of a calculation at gridsize Ax and viscous length L in the vicinity of the shock as 

V^Ax(x, t- L) = t-L)+r (^^, i; ^) , (C26) 

where the truncation error T depends on Ax through the ratio L/ /S.x, as discussed above. Then 
the deviation function for calculations with the same value of L/A.X reduces to 

D^'^'\m) = i^'~'Hxs + L'~'^^,t;LW)+TU,t;^ 

V Ax 

- ^(2) (^^ + i(2)^, ^(D) _ r 1. (C27) 
= V(i)(x, + L«e,t;L«)-^(2)(^^+^(2)^^^.^(i))^ 

that is, the deviation function for calculations without truncation errors. 

We use this result to determine Ax — > deviation functions over a wide range of 
viscous lengths L, including values of L for which sufficiently accurate calculations (i.e., with 
small truncation errors) are computationally infeasible. In Figure 20 we show the scalar 
deviation measure ( C24| ) as a function of L, determined from a sequence of calculations at 



gridsizes Ax = 1/100,1/200,... ,1/3200, with L/Ax = 4. In the reinterpretation of these 
calculations as a sequence in increasing time (rather than increasing resolution), we have a 
sequence at times t = 1/32, 1/16, . . . , 1. Because the deviations result from a diffusive relaxation 
process, we naively expect their amplitudes to decrease exponentially with increasing time (or 
decreasing L). Clearly the deviations decrease, but not exponentially. (The behavior appears 
more similar to a power law, but it does not approach one asymptotically. In any case, we have 
no reason to expect such behavior.) 



Morduchow and Paullay ( |1971| ) and Morduchow and Valentino ( 1979| ), working with a 



one-dimensional ideal gas with a Prandtl number of 3/4, performed an eigenvalue analysis to 
investigate shock stability. They found a continuum of decay rates, including values arbitrarily 
close to zero. If our system shares this property, then even asymptotically we should not expect to 
observe exponential decay. A combination of perturbation modes with a continuum of rates could 
imitate the observed behavior. 

Large-scale gradients (wall shock) 

Our third and last cause of deviations from exact scaling has a simple physical origin. 
Large-scale gradients in the flow have an effect on the structure of the transition. Since transitions 
of different widths sample these gradients to different extents, observable deviations result. To 
illustrate, we consider a physical problem in which the shock itself is not in steady state, the 
wall-shock problem. Although as t^oo the shock approaches steady state, at any finite time the 
shock generates weak gradients in the downstream region. 
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Figure 21 shows deviation functions D[p\{^) for the wall-shock problem. The two curves 
represent comparisons of, first (solid line), L^^) = 1/128 and L^^) = 1/256 (using L/Ax = 8), 
and second (dashed line), L^^^ = 1/256 and L^^) = 1/512 (using L/Ax = 4). (The sense of the 
deviation shows that the transition is steepening; this agrees with the smoothed initial conditions, 
which [by design] are less steep than the steady-state transition.) Note the effect of the large-scale 
gradients downstream from the shock — unlike in the steady-state, shock-tube cases, D{^) tends 
toward non-zero values for negative values of ^ = (x — Xs)/L. The slope (with respect to ^) with 
which the curve turns away from zero on the left side should be proportional to the scaling factor, 
i.e., L, as Figure 21 verifies. 



C.4. Extrapolation difficulties 



Our investigations into extrapolation procedures suggest difficulties in devising an automated 
extrapolation scheme that ensures a predetermined order of accuracy. Extrapolating sequences 
with the artificial-viscous length L fixed appears feasible, but the improvements arc small. Direct 
pointwise extrapolation of sequences in which the viscous length L is chosen proportional to the 
gridsize Ax yields improvements away from shocks, but at the cost of reasonable behavior at 
shocks. Finally, our studies of the scaling behavior of shock transitions suggests that computed 
transitions in general deviate from exact scaling behavior in ways not easily controlled; exploiting 
the scaling behavior to extrapolate to zero-width shock transitions may incur errors of unknown 
magnitude. 



REFERENCES 



Finn, L. S. and Hawley, J. F. 1989, unpublished report. 

Hawley, J. F., Smarr, L. L. and Wilson, J. R. 1984, ApJS, 55, 211. 

Landau, L. D. and Lifshitz, E. M. 1987, Fluid Mechanics, second edition, translated by J. B. Sykes 
and W. H. Reid (Oxford: Pergamon). 

Lax, R and Wendroff, B. 1960, Comm. Pure Appl. Math., 13, 217. 

Morduchow, M. and Paullay, A. J. 1971, Rhys. Fluids, 14, 323. 

Morduchow, M. and Valentino, J. V. 1979, J. Appl. Mech., 46, 505. 
Norman, M. L., Wilson, J. B. and Barton, R. T. 1980, ApJ, 239, 968. 

Norman, M. L. and Winkler, K.-H. A. 1986, 2-D Eulerian Hydrodynamics with Fluid Interfaces, 
Self-Gravity and Rotation, in Astrophysical Radiation Hydrodynamics, edited by M. L. 
Norman and K.-H. A. Winkler (Dordrecht: Reidel), p. 187. 

Richardson, L. F. and Gaunt, J. A. 1927, Phil. Trans. Roy. Soc. Lond. A, 226, 299. 

Richtmyer, R. D. and Morton, K. W. 1967, Difference Methods for Initial- Value Problems, second 
edition (New York: Interscience) . 

Sod, G. A. 1978, J. Comput. Phys., 27, 1. 

Sod, G. A. 1985, Numerical Methods in Fluid Dynamics: Initial and Initial Boundary- Value 
Problems (Cambridge: Cambridge University Press). 

Stone, J. M. and Norman, M. L. 1992, ApJS, 80, 753. 

van Leer, B. 1977, J. Comput. Phys., 23, 276. 



This preprint was prepared with the AAS lAT^^X macros v3.0. 



- 28 - 



Figure captions 

1. Difference between total energy in an internal-energy computation and exact total energy, 
as a function of gridsize Ax in the steepening-wave problem of section 3. Artificial viscosity 
vanishes with decreasing gridsize. Each line corresponds to a different time (from bottom to 
top, t = 4, 8, 12, . . . , 40). When t > 24, the error does not tend toward zero. 

2. Cauchy errors ( pl| ) for momentum density in the steepening-wave problem. The two lines 
correspond to t = 20 and t = 40. 

3. Difference between total energy in an internal-energy computation and exact total energy, as 
a function of gridsize Ax in the steepening-wave problem. The viscous lengths Li and L2 are 
fixed. Each line corresponds to a different time {t = 4, 8, 12, ... , 40). 

4. Cauchy errors for the steepening-wave problem. The viscous lengths Li and L2 are fixed. The 
two lines correspond to t = 20 and t = 40. 

5. The region near the shock front in the shock-tube problem calculated by both internal-energy 
(solid curve) and total-energy schemes (dashed curve) (at t = 1, using Ax = 1/6400 

and Li = L2 = 3/12800) as well as the exact solution (dotted curve). 

6. Difference between total energy in an internal-energy computation of the shock tube and the 
exact total energy, as a function of gridsize Ax. Increased resolution does not yield a better 
result. 

7. Momentum Cauchy errors for both internal-energy and total-energy calculations of the shock 
tube. Their behaviors are almost indistinguishable, although the internal-energy calculation is 
wrong. 

8. Momentum Cauchy errors for a small-amplitude sound wave. The quadratic viscous length L2 
is zero. The linear viscous length Li is chosen to be 4Ax. The short line segment indicates a 
slope corresponding to Ax dependence. 

9. The numerical convergence rate /3 for a small-amplitude sound wave. The viscous lengths 
are Li = 4Ax and L2 = 0. 

10. Cauchy errors for the analytic solution of the heat equation as diffusion coefficient i/— >0. As 
discussed in the text, the linearized equations of motion in the presence of linear artificial 
viscosity may take this form. 

11. The numerical convergence rate /5 associated with solutions of the heat equation as v^O. 

12. Density as a function of x for the wall-shock problem of section 5, at t = 2, computed 
with Li = L2 = 1/128 at Ax = 1/256. 
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13. Momentum Cauchy errors for the wall-shock problem. Each solid line describes a sequence of 
calculations in which the viscous length Li = L2 = L is a constant multiple of the gridsize; 
five lines represent the values L/Ax = 1,2,4,8, 16. The short line segment indicates a slope 
corresponding to Ax dependence. The dashed and dotted lines isolate the contributions from 
the shock-transition region {x < 7) and the rest of the flow (x > 7), respectively. 

14. Momentum Cauchy errors for the wall-shock problem. Each line describes a sequence of 
calculations in which the viscous lengths Li = L2 = L are held constant (from left to right, 
L = 1/256,1/128, ... ,1/8). 

15. Convergence of the momentum Cauchy errors for the wall-shock problem, considered within 
sequences in which the viscous lengths Li and L2 are constant (horizontal lines) and 
proportional to the gridsize Ax (diagonals). Numerical convergence rates /3, shown along these 
sequences, measure the approach to the expected convergence rates 0{Ax) and 0((Ax)^). 

16. Momentum Cauchy errors for extrapolations in the wall-shock problem, with L = 1/8. The 
extrapolations use calculations from Ax = 1/16 down to the value of the gridsize shown on the 
horizontal axis. 

17. The result of a pointwise extrapolation of three calculations of the wall-shock problem, at 
gridsizes Ax = 1/64, 1/128, 1/256, with L/Ax = 2. 

18. The 0{{Axf) truncation error (estimated as discussed in text) as a function of ^ = (x — Xs)/L 
in the shock-tube problem at t = 1. The estimates are for four calculations, all with L/Ax = 4, 
with gridsizes Ax = 1/100, 1/200, 1/400, 1/800. As discussed in the text, these calculations 
may also be interpreted as a sequence in time {t = 1/8, 1/4, 1/2, 1) at Ax = 1/800. 

19. Deviation from exact shock-transition scaling of density ( P23| ) as a function of ^ = (x — Xs)/L 
in the shock-tube problem. The solid line compares L = 1/25 and L = 1/50 calculations; 
the 0((Ax)^ truncation error is negligible. The dashed line compares L = 1/25, Ax = 1/100 
and L = 1/50, Ax = 1/200 calculations; the 0((Ax)^) truncation error for either is larger in 
magnitude than the function plotted here. 



20. Deviation from exact shock-transition scaling of density ( C24| ) as a function of viscous length L. 
The 0((Ax)^) truncation error has been removed, as described in the text, leaving the effects 
of relaxation toward exact scaling. 



21. Deviation from exact shock-transition scaling of density ( p23 ) as a function of = (x — Xs)/L 



for the wall-shock problem. The solid curve represents a comparison of L = 1/128 
with L = 1/256; the dashed curve, of L = 1/256 with L = 1/512. 
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